The following benchmark is of a 1122 ODEs with 24388 terms that describe a stiff chemical reaction network.
using ReactionNetworkImporters, OrdinaryDiffEq, DiffEqBiological, Sundials, Plots, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq, LSODA, TimerOutputs gr() prnbng = loadrxnetwork(BNGNetwork(), "BNGRepressilator", joinpath(dirname(pathof(ReactionNetworkImporters)),"..","data","bcr","bcr.net"))
Parsing parameters...done Adding parameters...done Parsing species...done Adding species...done Parsing and adding reactions...done Parsing groups...done
rn = deepcopy(prnbng.rn) addodes!(rn; build_jac=false, build_symfuncs=false, build_paramjac=false) tf = 100000.0 oprob = ODEProblem(rn, prnbng.u₀, (0.,tf), prnbng.p); densejac_rn = deepcopy(prnbng.rn) # zeroout_jac=true is needed to keep the generated expressions from being too big for the compiler addodes!(densejac_rn; build_jac=true, zeroout_jac = true, sparse_jac = false, build_symfuncs=false, build_paramjac=false) densejacprob = ODEProblem(densejac_rn, prnbng.u₀, (0.,tf), prnbng.p); sparsejac_rn = deepcopy(prnbng.rn) addodes!(sparsejac_rn; build_jac=true, sparse_jac = true, build_symfuncs=false, build_paramjac=false) sparsejacprob = ODEProblem(sparsejac_rn, prnbng.u₀, (0.,tf), prnbng.p);
@show numspecies(rn) # Number of ODEs
numspecies(rn) = 1122
@show numreactions(rn) # Apprx. number of terms in the ODE
numreactions(rn) = 24388
@show numparams(rn) # Number of Parameters
numparams(rn) = 128 128
As compiling the ODE derivative functions has in the past taken longer than running a simulation, we first force compilation by evaluating these functions one time.
const to = TimerOutput() u₀ = prnbng.u₀ u = copy(u₀); du = similar(u); p = prnbng.p @timeit to "ODERHS Eval1" rn.f(du,u,p,0.) @timeit to "ODERHS Eval2" rn.f(du,u,p,0.) sparsejac_rn.f(du,u,p,0.) J = zeros(length(u),length(u)) @timeit to "DenseJac Eval1" densejac_rn.jac(J,u,p,0.) @timeit to "DenseJac Eval2" densejac_rn.jac(J,u,p,0.) Js = similar(sparsejac_rn.odefun.jac_prototype) @timeit to "SparseJac Eval1" sparsejac_rn.jac(Js,u,p,0.) @timeit to "SparseJac Eval2" sparsejac_rn.jac(Js,u,p,0.) show(to)
──────────────────────────────────────────────────────────────────────────
Time Allocations
────────────────────── ───────────────────────
Tot / % measured: 127s / 84.7% 17.2GiB / 80.2%
Section ncalls time %tot avg alloc %tot avg
──────────────────────────────────────────────────────────────────────────
DenseJac Eval1 1 52.5s 49.0% 52.5s 5.29GiB 38.3% 5.29GiB
SparseJac Eval1 1 35.6s 33.3% 35.6s 5.23GiB 37.8% 5.23GiB
ODERHS Eval1 1 19.0s 17.7% 19.0s 3.31GiB 23.9% 3.31GiB
DenseJac Eval2 1 573μs 0.00% 573μs 32.0B 0.00% 32.0B
ODERHS Eval2 1 84.4μs 0.00% 84.4μs 32.0B 0.00% 32.0B
SparseJac Eval2 1 49.7μs 0.00% 49.7μs 32.0B 0.00% 32.0B
──────────────────────────────────────────────────────────────────────────
sol = solve(oprob, CVODE_BDF(), saveat=tf/1000., reltol=1e-5, abstol=1e-5) plot(sol,legend=false, fmt=:png)